Step 0: Prepare covariates and input files

IDs <- fread("~/genotype_qc/TERRE_QC/all_imputed_r2_30_rsid_hard_call.fam")[, .(FID = V1, IID = V2)]
prsice_cov <- fread("prsice_cov_and_status_mvalues.txt")
prsice_cov <- prsice_cov[match(IDs$IID,prsice_cov$IID)]
all(IDs$IID ==prsice_cov$IID)
[1] NA
covariate <- cbind(IDs, prsice_cov[,-c(16,17,18)])

PD <- cbind(IDs,prsice_cov[,16])

head(covariate)
head(PD)
fwrite(na.omit(PD), "TERRE.pheno", sep = "\t")
fwrite(na.omit(covariate), "TERRE.covariate", sep = "\t")
covariate

Step 1: Run PRSice-2 on Nalls et al 2019 Sumstats

Rscript /home1/NEURO/casazza/PRSice.R \
    --prsice /home1/NEURO/casazza/PRSice_linux\
    --base /home1/NEURO/casazza/nalls_PD.QC.gz\
    --base-info INFO:0.8 \
    --base-maf MAF:0.01 \
    --cov TERRE.covariate \
    --binary-target T\
    --beta  \
    --ld /home1/NEURO/casazza/1000G_plink/EUR_phase3  \
    --out TERRE_PRSice \
    -q 5\
    --all-score\
    --pheno TERRE.pheno \
    --snp SNP \
    --stat b \
    --pvalue p\
    --target /home1/NEURO/casazza/genotype_qc/TERRE_QC/all_imputed_r2_30_rsid_hard_call \
    --thread 32

Step 2: Evaluate output

include_graphics("prsice_images/TERRE_PRSice_BARPLOT_2022-06-30.png")

include_graphics("prsice_images/TERRE_PRSice_HIGH-RES_PLOT_2022-06-30.png")

include_graphics("prsice_images/TERRE_PRSice_QUANTILES_PLOT_2022-06-30.png")


include_graphics("prsice_images/TERRE_PRSice_nalls_male_BARPLOT_2022-10-04.png")

include_graphics("prsice_images/TERRE_PRSice_nalls_male_HIGH-RES_PLOT_2022-10-04.png")

include_graphics("prsice_images/TERRE_PRSice_nalls_male_QUANTILES_PLOT_2022-10-04.png")


include_graphics("prsice_images/TERRE_PRSice_nalls_female_BARPLOT_2022-10-04.png")

include_graphics("prsice_images/TERRE_PRSice_nalls_female_HIGH-RES_PLOT_2022-10-04.png")

include_graphics("prsice_images/TERRE_PRSice_nalls_female_QUANTILES_PLOT_2022-10-04.png")

Plotting PRSice Data on my own

library(ggnewscale)
prsice_male_meta <- fread("prsice_nalls_male_data/TERRE_PRSice_nalls_male.prsice")
ggplot(prsice_male_meta[Threshold <= 0.5], aes(Threshold, Num_SNP, color = -log10(P))) +
  geom_point() +
  scale_y_continuous(breaks = c(seq(0, 1e5, 2.5e4), seq(2e5, 6e5, 1e5))) +
  theme_minimal()


prsice_female_meta <- fread("prsice_nalls_female_data/TERRE_PRSice_nalls_female.prsice")
ggplot(prsice_female_meta[Threshold <= 0.5], aes(Threshold, Num_SNP, color = -log10(P))) +
  geom_point() +
  scale_y_continuous(breaks = c(seq(0, 1e5, 2.5e4), seq(2e5, 6e5, 1e5))) +
  theme_minimal()


prsice_meta <- fread("prsice_data/TERRE_PRSice.prsice")
ggplot(mapping = aes(Threshold, R2, color = -log10(P))) +
  geom_point(data = prsice_female_meta, size = 1) +
  scale_color_gradient(low = "lightpink4", high = "lightpink") +
  labs(color = bquote("Female log"["10"] ~ "(P)")) +
  new_scale_color() +
  geom_point(data = prsice_male_meta, size = 1, aes(color = -log10(P))) +
  scale_color_gradient(low = "lightblue4", high = "lightblue") +
  labs(y = bquote("R"^2), x = "GWAS P-Value Threshold", color = bquote("Male -log"["10"] ~ "(P)")) +
  theme_minimal() +
  theme(legend.position = "top")

ggplot(mapping = aes(Threshold, R2, color = -log10(P))) +
  geom_point(data = prsice_female_meta, size = 1) +
  scale_color_gradient(low = "lightpink4", high = "lightpink",guide = guide_colorbar(order=3)) +
  labs(color = bquote("Female -log"["10"] ~ "(P)")) +
  new_scale_color() +
  geom_point(data = prsice_male_meta, size = 1, aes(color = -log10(P))) +
  scale_color_gradient(low = "lightblue4", high = "lightblue",guide = guide_colorbar(order=2)) +
  labs(color = bquote("Male -log"["10"] ~ "(P)")) +
  new_scale_color() +
  geom_point(data = prsice_meta, size = 1, aes(color = -log10(P))) +
  scale_color_gradient(low = "gray40", high = "gray80",guide = guide_colorbar(order=1)) +
  labs(y = bquote("R"^2), x = "GWAS P-Value Threshold", color = bquote("Cross-sex -log"["10"] ~ "(P)")) +
  theme_minimal() +
  theme(legend.position = "top",legend.title = element_text(size=7))

Step 3 run linear model at different thresholds for SNP inclusion

This is updated to recently processed DNAm Data 20 outliers in DNAm removed:

load("/home1/NEURO/SHARE_DECIPHER/processed_DNAm_data/2022/TERRE_processed_2022/1-TERRE_RG_filtered.RData") #PD_RG_filtered
# Assign genotyping ID to data
original_covars <- fread("/home1/NEURO/SHARE_DECIPHER/terre_meta_master.csv")[, .(patient, IID = gsub("_PAE.*", "", IID))]
betas_combat <- minfi::getBeta(PD_RG_filtered)
colnames(betas_combat) <- original_covars$IID[match(colnames(betas_combat), original_covars$patient)]
betas_combat <- betas_combat[, colnames(betas_combat) %in% covariate$IID]

Let’s check how the data looks for the first 5 subjects:

ggplot(betas_combat[, 1:5] %>% as.data.table(keep.rownames = T) %>% melt(id.vars = "rn", value.name = "betas", variable.name = "subject"), aes(betas, color = subject))+
  geom_density()

Match DNA, PRS, and metadata @TODO fix this after update of PRS etc.

# prsice_best=0.00030005
prsice_all <- fread("prsice_data/TERRE_PRSice.all_score")[match(colnames(betas_combat), IID,nomatch=0), .(FID, IID, `Pt_5e-08`, `Pt_5.005e-05`, `Pt_0.00010005`, `Pt_0.00100005`, `Pt_0.00030005`,`Pt_0.0101501`, `Pt_0.1`, `Pt_0.2`, `Pt_0.3`, `Pt_0.4`, `Pt_0.5`, `Pt_1`)]
covariate <- covariate[match(colnames(betas_combat), IID)]
all(covariate$IID == colnames(betas_combat))
[1] TRUE
all(covariate$IID == prsice_all$IID)
[1] TRUE
covariate_male <- covariate[sex == 1] %>% select(-sex)
betas_male <- betas_combat[, covariate_male$IID]
# best = 0.0335001
prsice_male_all <-fread("prsice_nalls_male_data/TERRE_PRSice_nalls_male.all_score")[match(colnames(betas_male), IID,nomatch=0), .(FID, IID, `Pt_5e-08`, `Pt_5.005e-05`, `Pt_0.00010005`, `Pt_0.00100005`, `Pt_0.0101501`, `Pt_0.0335001`, `Pt_0.1`, `Pt_0.2`, `Pt_0.3`, `Pt_0.4`, `Pt_0.5`, `Pt_1`)]

covariate_male <- covariate_male[match(colnames(betas_male), IID,nomatch=0)]
all(covariate_male$IID == colnames(betas_male))
[1] TRUE
all(covariate_male$IID == prsice_male_all$IID)
[1] TRUE
covariate_female <- covariate[sex == 0] %>% select(-sex) #%>% filter(IID %in% prsice_female_best$IID)
betas_female <- betas_combat[, covariate_female$IID]
#best = 0.00035005

prsice_female_all <- fread("prsice_nalls_female_data/TERRE_PRSice_nalls_female.all_score")[match(colnames(betas_female), IID,nomatch=0), .(FID, IID, `Pt_5e-08`, `Pt_5.005e-05`, `Pt_0.00010005`, `Pt_0.00035005`, `Pt_0.00100005`, `Pt_0.0101501`, `Pt_0.1`, `Pt_0.2`, `Pt_0.3`, `Pt_0.4`, `Pt_0.5`, `Pt_1`)]
covariate_female <- covariate_female[match(colnames(betas_female), IID,nomatch=0)]
all(covariate_female$IID == colnames(betas_female))
[1] TRUE
all(covariate_female$IID == prsice_female_all$IID)
[1] TRUE

Run limma

mvalues <- lumi::beta2m(betas_combat)
prs_mat <- prsice_all[, -c(1, 2)]
cov_mat <- covariate[, -c(1, 2)]

mvalues_male <- lumi::beta2m(betas_male)
prs_mat_male <- prsice_male_all[, -c(1, 2)]
cov_mat_male <- covariate_male[, -c(1, 2)]

mvalues_female <- lumi::beta2m(betas_female)
prs_mat_female <- prsice_female_all[, -c(1, 2)]
cov_mat_female <- covariate_female[, -c(1, 2)]
registerDoParallel(ncol(prs_mat) / 4)
hits <- foreach(prs_thresh = colnames(prs_mat)) %dopar% {
  design_prs <- model.matrix(~., data = cbind(prs_mat[, ..prs_thresh], cov_mat))
  prs_fit <- lmFit(mvalues, design_prs)
  prs_fit <- eBayes(prs_fit)
  topTable(prs_fit, coef = 2, adjust.method = "bonf", p.value = 0.05, number = Inf, genelist = rownames(mvalues))
}
names(hits) <- colnames(prs_mat)
hits_by_thresh_bonf <- rbindlist(hits, idcol = "threshold", fill = TRUE)

registerDoParallel(ncol(prs_mat_male) / 4)
hits_male <- foreach(prs_thresh = colnames(prs_mat_male)) %dopar% {
  design_prs_male <- model.matrix(~., data = cbind(prs_mat_male[, ..prs_thresh], cov_mat_male))
  prs_fit_male <- lmFit(mvalues_male, design_prs_male)
  prs_fit_male <- eBayes(prs_fit_male)
  topTable(prs_fit_male, coef = 2, adjust.method = "bonf", p.value = 0.05, number = Inf, genelist = rownames(mvalues_male))
}
names(hits_male) <- colnames(prs_mat_male)
hits_by_thresh_bonf_male <- rbindlist(hits_male, idcol = "threshold", fill = TRUE)

registerDoParallel(ncol(prs_mat_female) / 4)
hits_female <- foreach(prs_thresh = colnames(prs_mat_female)) %dopar% {
  design_prs_female <- model.matrix(~., data = cbind(prs_mat_female[, ..prs_thresh], cov_mat_female))
  prs_fit_female <- lmFit(mvalues_female, design_prs_female)
  prs_fit_female <- eBayes(prs_fit_female)
  topTable(prs_fit_female, coef = 2, adjust.method = "bonf", p.value = 0.05, number = Inf, genelist = rownames(mvalues_female))
}

Plotting EWAS vs Threshold Experiment by Sex

ggplot(to_plot, aes(threshold, hits, fill = Sex, label = hits)) +
  geom_text(position = plot_pos, vjust = -0.25) +
  geom_col(position = plot_pos) +
  labs(x = "GWAS P Value Threshold", y = "EWAS Hits") +
  scale_fill_manual(values = c("grey80", "lightblue", "lightpink")) +
  theme_minimal()
hits_by_thresh_bonf[, .(hits = .N, Sex = "Cross-sex"), by = threshold]

hits_by_thresh_bonf_male[, .(hits = .N), by = threshold]
hits_by_thresh_bonf_female[, .(hits = .N), by = threshold]
display_venn <- function(x, ...) {
  library(VennDiagram)
  grid.newpage()
  venn_object <- venn.diagram(x, filename = NULL, ...)
  grid.draw(venn_object)
}

display_venn(list(`Cross-sex` = hits_by_thresh_bonf[threshold == "Pt_5e-08"]$ID, Male = hits_by_thresh_bonf_male[threshold == "Pt_5e-08"]$ID, Female = hits_by_thresh_bonf_female[threshold == "Pt_5e-08"]$ID), fill = c("gray80", "lightblue", "lightpink"))
Loading required package: grid
Loading required package: futile.logger

male_ids <- hits_by_thresh_bonf_male[threshold == "Pt_5e-08"]$ID
cross_ids <- hits_by_thresh_bonf[threshold == "Pt_5e-08"]$ID
male_ids[!male_ids %in% cross_ids]
character(0)
get_full_fit <- function(prs_mat,cov_mat,mvalues){
  top_design_prs <- model.matrix(~., data = cbind(prs_mat[, `Pt_5e-08`], cov_mat))
  top_prs_fit <- lmFit(mvalues, top_design_prs)
  top_prs_fit <- eBayes(top_prs_fit)
  top_prs_hits <- topTable(top_prs_fit, coef = 2, adjust.method = "bonf", number = Inf, genelist = rownames(mvalues))
}
top_prs_hits <- get_full_fit(prs_mat,cov_mat,mvalues)
top_male_prs_hits <- get_full_fit(prs_mat_male, cov_mat_male, mvalues_male)
top_female_prs_hits <- get_full_fit(prs_mat_female, cov_mat_female, mvalues_female)
save(list=c("top_prs_hits","top_male_prs_hits","top_female_prs_hits"),file="prs_nalls_cross_w_sex_stratified.RData")
load("prs_nalls_cross_w_sex_stratified.RData")
manifest <- IlluminaHumanMethylationEPICanno.ilm10b4.hg19::Other %>%
  as.data.frame() %>%
  rownames_to_column(var = "name")
prs_annot <- data.table(top_prs_hits)[manifest, gene := gsub(";.*", "", UCSC_RefGene_Name), on = c(ID = "name")]
prs_annot_male <- data.table(top_male_prs_hits)[manifest, gene := gsub(";.*", "", UCSC_RefGene_Name), on = c(ID = "name")]
prs_annot_female<- data.table(top_female_prs_hits)[manifest, gene := gsub(";.*", "", UCSC_RefGene_Name), on = c(ID = "name")]
plot_prs_hits <- function(prs_annot,label_color){
  ggplot(prs_annot, aes(logFC, -log10(P.Value))) +
    geom_point() +
    geom_point(
      data = subset(prs_annot, adj.P.Val < 0.05 & abs(logFC) > 0.03),
      color = label_color,
      mapping = aes(logFC, -log10(P.Value))
    ) +
    geom_hline(
      linetype = "dashed",
      yintercept = min(-log10(prs_annot$P.Value[prs_annot$adj.P.Val < 0.05]))
    ) +
    geom_vline(linetype = "dashed", xintercept = 0.03) +
    geom_vline(linetype = "dashed", xintercept = -0.03) +
    geom_text_repel(
      data = prs_annot %>% filter(abs(logFC) > 0.03 & adj.P.Val < 0.05),
      color = "dodgerblue",
      mapping = aes(logFC, -log10(P.Value), label = ifelse(gene != "", gene, ID)),
      size = 3,
      max.overlaps = 20
    ) +
    labs(y = bquote("log"[10] ~ "(P)"), x = quote(Delta ~ "M" ~ Methylation)) +
    theme_minimal()
}

plot_prs_hits(prs_annot,"gray40")

plot_prs_hits(prs_annot_male,"lightblue")

plot_prs_hits(prs_annot_female,"pink")

go enrichment?

library(gprofiler2)
cross_sex <- unique(gsub(";.*","",manifest[manifest$name %in% hits_by_thresh_bonf[threshold == "Pt_5e-08"]$ID,]$UCSC_RefGene_Name))
background <- unique(gsub(";.*","",manifest$UCSC_RefGene_Name))
gost_res <- gost(query=cross_sex,custom_bg = background)
gostplot(gost_res)
gmirror(
  prs_annot_male[,.(SNP=ID,CHR=gsub("chr","",cpg_pos[match(ID,cpg_pos$name),]$chr),POS=cpg_pos[match(ID,cpg_pos$name),]$pos,pvalue=P.Value)],
  prs_annot_female[,.(SNP=ID,CHR=gsub("chr","",cpg_pos[match(ID,cpg_pos$name),]$chr),POS=cpg_pos[match(ID,cpg_pos$name),]$pos,pvalue=P.Value)],
  annotate_snp = c(male_annot_cpg,female_annot_cpg),
  tline =  max(prs_annot_male[adj.P.Val<0.25]$P.Value),
  bline =  max(prs_annot_female[adj.P.Val<0.25]$P.Value),
  highlight_p = c(max(prs_annot_male[adj.P.Val<0.25]$P.Value),max(prs_annot_female[adj.P.Val<0.25]$P.Value)),
  toptitle="Male",
  bottomtitle ="Female",
  highlighter="green",
  background="white"
)
[1] "Saving plot to gmirror.png"
TableGrob (2 x 1) "arrange": 2 grobs

Manhattan plot vs GWAS manhattan plot

library(qqman)

copy_annot <- prs_annot[cpg_pos, on = c(ID = "name")] %>%
  mutate(chr = as.numeric(recode(gsub("chr", "", chr),X="23",Y="24")))
to_plot <- copy_annot[, .(SNP = gene, CHR = chr, BP = pos, P = P.Value, FDR = adj.P.Val)][!is.na(P)]
manhattan(to_plot,
  annotatePval = max(to_plot[FDR < 0.05]$P),
  annotateTop = TRUE,
  genomewideline = min(-log10(to_plot[FDR < 0.05]$P)),
  suggestiveline = FALSE
)
manhattan(to_plot[CHR == 17],
  annotatePval = max(to_plot[FDR < 0.05]$P),
  annotateTop = FALSE,
  genomewideline = min(-log10(to_plot[FDR < 0.05]$P)),
  suggestiveline = FALSE
)
qq(to_plot$P)

DMRs

library(DMRcate)
Loading required package: minfi
Loading required package: BiocGenerics

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, clusterExport, clusterMap, parApply, parCapply,
    parLapply, parLapplyLB, parRapply, parSapply, parSapplyLB

The following object is masked from ‘package:limma’:

    plotMA

The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated, eval, evalq,
    Filter, Find, get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget, order, paste, pmax,
    pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply, union,
    unique, unsplit, which, which.max, which.min

Loading required package: GenomicRanges
Loading required package: stats4
Loading required package: S4Vectors

Attaching package: ‘S4Vectors’

The following objects are masked from ‘package:data.table’:

    first, second

The following objects are masked from ‘package:dplyr’:

    first, rename

The following object is masked from ‘package:tidyr’:

    expand

The following object is masked from ‘package:base’:

    expand.grid

Loading required package: IRanges

Attaching package: ‘IRanges’

The following object is masked from ‘package:data.table’:

    shift

The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice

The following object is masked from ‘package:purrr’:

    reduce

Loading required package: GenomeInfoDb
Loading required package: SummarizedExperiment
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: DelayedArray
Loading required package: matrixStats

Attaching package: ‘matrixStats’

The following objects are masked from ‘package:Biobase’:

    anyMissing, rowMedians

The following object is masked from ‘package:dplyr’:

    count

Loading required package: BiocParallel

Attaching package: ‘DelayedArray’

The following objects are masked from ‘package:matrixStats’:

    colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges

The following object is masked from ‘package:purrr’:

    simplify

The following objects are masked from ‘package:base’:

    aperm, apply, rowsum

Loading required package: Biostrings
Loading required package: XVector

Attaching package: ‘XVector’

The following object is masked from ‘package:purrr’:

    compact


Attaching package: ‘Biostrings’

The following object is masked from ‘package:base’:

    strsplit

Loading required package: bumphunter
Loading required package: locfit
locfit 1.5-9.4   2020-03-24

Attaching package: ‘locfit’

The following object is masked from ‘package:purrr’:

    none

replacing previous import 'minfi::getMeth' by 'bsseq::getMeth' when loading 'DMRcate'
S4_to_dataframe <- function(s4obj) {
  nms <- slotNames(s4obj)
  lst <- lapply(nms, function(nm) slot(s4obj, nm))
  as.data.frame(setNames(lst, nms))
}
run_dmrcate <- function(prs_mat,cov_mat,mvalues){
  design_prs <- model.matrix(~., data = cbind(prs_mat[, `Pt_5e-08`], cov_mat))
  prs_annotated <- cpg.annotate(datatype = "array", object = mvalues, analysis.type = "differential", design = design_prs, coef = 2, what = "M", arraytype = "EPIC", fdr = 0.05)
  prs_dmr_res <- dmrcate(prs_annotated, lambda = 1000, C = 2)
  return(S4_to_dataframe(prs_dmr_res))
}
dmr_cross <- run_dmrcate(prs_mat, cov_mat, mvalues)
Loading required package: IlluminaHumanMethylationEPICanno.ilm10b4.hg19
Your contrast returned 85 individually significant probes; a small but real effect. Consider manually setting the value of pcutoff to return more DMRs, but be warned that doing this increases the risk of Type I errors.
Fitting chr1...
Fitting chr2...
Fitting chr3...
Fitting chr4...
Fitting chr5...
Fitting chr6...
Fitting chr7...
Fitting chr8...
Fitting chr9...
Fitting chr10...
Fitting chr11...
Fitting chr12...
Fitting chr13...
Fitting chr14...
Fitting chr15...
Fitting chr16...
Fitting chr17...
Fitting chr18...
Fitting chr19...
Fitting chr20...
Fitting chr21...
Fitting chr22...
Fitting chrX...
Fitting chrY...
Demarcating regions...
Done!
dmr_males <- run_dmrcate(prs_mat_male,cov_mat_male,mvalues_male)
Your contrast returned 59 individually significant probes; a small but real effect. Consider manually setting the value of pcutoff to return more DMRs, but be warned that doing this increases the risk of Type I errors.
Fitting chr1...
Fitting chr2...
Fitting chr3...
Fitting chr4...
Fitting chr5...
Fitting chr6...
Fitting chr7...
Fitting chr8...
Fitting chr9...
Fitting chr10...
Fitting chr11...
Fitting chr12...
Fitting chr13...
Fitting chr14...
Fitting chr15...
Fitting chr16...
Fitting chr17...
Fitting chr18...
Fitting chr19...
Fitting chr20...
Fitting chr21...
Fitting chr22...
Fitting chrX...
Fitting chrY...
Demarcating regions...
Done!
dmr_female <- run_dmrcate(prs_mat_female,cov_mat_female,mvalues_female)
Your contrast returned 41 individually significant probes; a small but real effect. Consider manually setting the value of pcutoff to return more DMRs, but be warned that doing this increases the risk of Type I errors.
Fitting chr1...
Fitting chr2...
Fitting chr3...
Fitting chr4...
Fitting chr5...
Fitting chr6...
Fitting chr7...
Fitting chr8...
Fitting chr9...
Fitting chr10...
Fitting chr11...
Fitting chr12...
Fitting chr13...
Fitting chr14...
Fitting chr15...
Fitting chr16...
Fitting chr17...
Fitting chr18...
Fitting chr19...
Fitting chr20...
Fitting chr21...
Fitting chr22...
Fitting chrX...
Fitting chrY...
Demarcating regions...
Done!
annotation <- minfi::getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)

Plotting all DMRs

get_dmr_effects <- function(dmr, limma_res,mvalues) {
  dmr_coord <- str_match_all(as.character(dmr), "(chr.*):([0-9]*)-([0-9]*)")[[1]]
  cpgs <- as.data.table(annotation[annotation$chr == dmr_coord[2] & annotation$pos >= as.numeric(dmr_coord[3]) & annotation$pos <= as.numeric(dmr_coord[4]), ])
  res <- limma_res[cpgs,on=c("ID"="Name"),nomatch=0 ]
  dmr_1 <- as.data.frame(res[!is.na(res$logFC), ])
  dmr_methy <- reshape2::melt(lumi::m2beta(mvalues[dmr_1$ID, ]), stringsAsFactors = FALSE)
  to_plot <- merge(dmr_1, dmr_methy, by.x = "ID", by.y = "Var1")
  to_plot$DMR <- dmr
  to_plot
}
get_dmr_res <- function(dmr, limma_res) {
  dmr_coord <- str_match_all(as.character(dmr), "(chr.*):([0-9]*)-([0-9]*)")[[1]]
  cpgs <- annotation[annotation$chr == dmr_coord[2] & annotation$pos >= as.numeric(dmr_coord[3]) & annotation$pos <= as.numeric(dmr_coord[4]), ]
  res <- cbind(cpgs, limma_res[cpgs$Name, ])
  res$DMR <- dmr
  return(res)
}
plot_dmrs <- function(dmr_res,limma_res,prs_mat, mvalues, case_control) {
  to_plot <- rbindlist(lapply(dmr_res$coord, function(dmr) get_dmr_effects(dmr, limma_res,mvalues))) %>%
    mutate(SCORE1_AVG = scale(prs_mat[match(Var2,IDs$IID),`Pt_5e-08`]),PD = PD[match(Var2,IID)]$PD) %>%
    filter(!is.na(SCORE1_AVG))
  for(cur_dmr in unique(to_plot$DMR)){
    cur_plot <- to_plot[DMR == cur_dmr]
    p1 <- ggplot(cur_plot , aes(pos, value, color = SCORE1_AVG)) +
      geom_point() +
      scale_color_gradient(low = rev(case_control)[1], high = rev(case_control)[2]) +
      theme_minimal() +
      labs(title = unique(cur_plot$DMR),x="POS",y=bquote("Methylation"~beta),color = "Normalized PD PRS")+
      theme(axis.text.x = element_text(angle = 90))
    p2 <- ggplot(cur_plot %>% mutate(PD = ifelse(PD == 1, "CASE", "CONTROL")), aes(factor(pos), value, color = PD)) +
      geom_boxplot(position = position_dodge(0.75)) +
      scale_color_manual(values = case_control) +
      theme_minimal() +
      stat_summary(aes(group = PD), fun = mean, geom = "line") +
      labs(title = unique(cur_plot$DMR),x="POS",y=bquote("Methylation"~beta),color = "PD status") +
      theme(axis.text.x = element_text(angle = 90))
    print(p1)
    print(p2)
  }
}

plot_dmrs(dmr_cross, prs_annot, prs_mat,mvalues, rev(c("gray80", "gray40")))

plot_dmrs(dmr_males,prs_annot_male,prs_mat_male, mvalues_male,rev(c("light blue", "lightblue4")))

plot_dmrs(dmr_female, prs_annot_female, prs_mat_female, mvalues_female,rev(c("lightpink", "lightpink4")))

dmr_cross
dmr_males
dmr_female
save(list=c("dmr_cross","dmr_males","dmr_female"),file="prs_dmr_nalls.RData")
LS0tCnRpdGxlOiAiUFJTaWNlMiBEZXZlbG9wbWVudCBvZiByaXNrIHNjb3JlIgpvdXRwdXQ6CiAgaHRtbF9ub3RlYm9vazoKICAgIHRvYzogeWVzCi0tLQoKYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0V9CmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KGdncmVwZWwpCmxpYnJhcnkoZGF0YS50YWJsZSkKbGlicmFyeShrbml0cikKbGlicmFyeShsaW1tYSkKbGlicmFyeShmb3JlYWNoKQpsaWJyYXJ5KGRvUGFyYWxsZWwpClN5cy5zZXRsb2NhbGUoIkxDX01FU1NBR0VTIiwgImVuX1VTLnV0ZjgiKQprbml0cjo6b3B0c19jaHVuayRzZXQoZWNobyA9IFRSVUUpCmBgYAoKIyBTdGVwIDA6IFByZXBhcmUgY292YXJpYXRlcyBhbmQgaW5wdXQgZmlsZXMKYGBge3J9CklEcyA8LSBmcmVhZCgifi9nZW5vdHlwZV9xYy9URVJSRV9RQy9hbGxfaW1wdXRlZF9yMl8zMF9yc2lkX2hhcmRfY2FsbC5mYW0iKVssIC4oRklEID0gVjEsIElJRCA9IFYyKV0KcHJzaWNlX2NvdiA8LSBmcmVhZCgicHJzaWNlX2Nvdl9hbmRfc3RhdHVzX212YWx1ZXMudHh0IikKcHJzaWNlX2NvdiA8LSBwcnNpY2VfY292W21hdGNoKElEcyRJSUQscHJzaWNlX2NvdiRJSUQpXQphbGwoSURzJElJRCA9PXByc2ljZV9jb3YkSUlEKQpjb3ZhcmlhdGUgPC0gY2JpbmQoSURzLCBwcnNpY2VfY292WywtYygxNiwxNywxOCldKQoKUEQgPC0gY2JpbmQoSURzLHByc2ljZV9jb3ZbLDE2XSkKCmhlYWQoY292YXJpYXRlKQpoZWFkKFBEKQpmd3JpdGUobmEub21pdChQRCksICJURVJSRS5waGVubyIsIHNlcCA9ICJcdCIpCmZ3cml0ZShuYS5vbWl0KGNvdmFyaWF0ZSksICJURVJSRS5jb3ZhcmlhdGUiLCBzZXAgPSAiXHQiKQpgYGAKYGBge3J9CmNvdmFyaWF0ZQpgYGAKCgojIFN0ZXAgMTogUnVuIFBSU2ljZS0yIG9uIE5hbGxzIGV0IGFsIDIwMTkgU3Vtc3RhdHMKYGBge2Jhc2gsZXZhbD1GQUxTRX0KUnNjcmlwdCAvaG9tZTEvTkVVUk8vY2FzYXp6YS9QUlNpY2UuUiBcCiAgICAtLXByc2ljZSAvaG9tZTEvTkVVUk8vY2FzYXp6YS9QUlNpY2VfbGludXhcCiAgICAtLWJhc2UgL2hvbWUxL05FVVJPL2Nhc2F6emEvbmFsbHNfUEQuUUMuZ3pcCiAgICAtLWJhc2UtaW5mbyBJTkZPOjAuOCBcCiAgICAtLWJhc2UtbWFmIE1BRjowLjAxIFwKICAgIC0tY292IFRFUlJFLmNvdmFyaWF0ZSBcCiAgICAtLWJpbmFyeS10YXJnZXQgVFwKICAgIC0tYmV0YSAgXAogICAgLS1sZCAvaG9tZTEvTkVVUk8vY2FzYXp6YS8xMDAwR19wbGluay9FVVJfcGhhc2UzICBcCiAgICAtLW91dCBURVJSRV9QUlNpY2UgXAogICAgLXEgNVwKICAgIC0tYWxsLXNjb3JlXAogICAgLS1waGVubyBURVJSRS5waGVubyBcCiAgICAtLXNucCBTTlAgXAogICAgLS1zdGF0IGIgXAogICAgLS1wdmFsdWUgcFwKICAgIC0tdGFyZ2V0IC9ob21lMS9ORVVSTy9jYXNhenphL2dlbm90eXBlX3FjL1RFUlJFX1FDL2FsbF9pbXB1dGVkX3IyXzMwX3JzaWRfaGFyZF9jYWxsIFwKICAgIC0tdGhyZWFkIDMyCmBgYAojIFN0ZXAgMjogRXZhbHVhdGUgb3V0cHV0CmBgYHtyLCBvdXQud2lkdGg9IjQwMHB4In0KaW5jbHVkZV9ncmFwaGljcygicHJzaWNlX2ltYWdlcy9URVJSRV9QUlNpY2VfQkFSUExPVF8yMDIyLTA2LTMwLnBuZyIpCmluY2x1ZGVfZ3JhcGhpY3MoInByc2ljZV9pbWFnZXMvVEVSUkVfUFJTaWNlX0hJR0gtUkVTX1BMT1RfMjAyMi0wNi0zMC5wbmciKQppbmNsdWRlX2dyYXBoaWNzKCJwcnNpY2VfaW1hZ2VzL1RFUlJFX1BSU2ljZV9RVUFOVElMRVNfUExPVF8yMDIyLTA2LTMwLnBuZyIpCgppbmNsdWRlX2dyYXBoaWNzKCJwcnNpY2VfaW1hZ2VzL1RFUlJFX1BSU2ljZV9uYWxsc19tYWxlX0JBUlBMT1RfMjAyMi0xMC0wNC5wbmciKQppbmNsdWRlX2dyYXBoaWNzKCJwcnNpY2VfaW1hZ2VzL1RFUlJFX1BSU2ljZV9uYWxsc19tYWxlX0hJR0gtUkVTX1BMT1RfMjAyMi0xMC0wNC5wbmciKQppbmNsdWRlX2dyYXBoaWNzKCJwcnNpY2VfaW1hZ2VzL1RFUlJFX1BSU2ljZV9uYWxsc19tYWxlX1FVQU5USUxFU19QTE9UXzIwMjItMTAtMDQucG5nIikKCmluY2x1ZGVfZ3JhcGhpY3MoInByc2ljZV9pbWFnZXMvVEVSUkVfUFJTaWNlX25hbGxzX2ZlbWFsZV9CQVJQTE9UXzIwMjItMTAtMDQucG5nIikKaW5jbHVkZV9ncmFwaGljcygicHJzaWNlX2ltYWdlcy9URVJSRV9QUlNpY2VfbmFsbHNfZmVtYWxlX0hJR0gtUkVTX1BMT1RfMjAyMi0xMC0wNC5wbmciKQppbmNsdWRlX2dyYXBoaWNzKCJwcnNpY2VfaW1hZ2VzL1RFUlJFX1BSU2ljZV9uYWxsc19mZW1hbGVfUVVBTlRJTEVTX1BMT1RfMjAyMi0xMC0wNC5wbmciKQpgYGAKIyMgUGxvdHRpbmcgUFJTaWNlIERhdGEgb24gbXkgb3duCmBgYHtyfQpsaWJyYXJ5KGdnbmV3c2NhbGUpCnByc2ljZV9tYWxlX21ldGEgPC0gZnJlYWQoInByc2ljZV9uYWxsc19tYWxlX2RhdGEvVEVSUkVfUFJTaWNlX25hbGxzX21hbGUucHJzaWNlIikKZ2dwbG90KHByc2ljZV9tYWxlX21ldGFbVGhyZXNob2xkIDw9IDAuNV0sIGFlcyhUaHJlc2hvbGQsIE51bV9TTlAsIGNvbG9yID0gLWxvZzEwKFApKSkgKwogIGdlb21fcG9pbnQoKSArCiAgc2NhbGVfeV9jb250aW51b3VzKGJyZWFrcyA9IGMoc2VxKDAsIDFlNSwgMi41ZTQpLCBzZXEoMmU1LCA2ZTUsIDFlNSkpKSArCiAgdGhlbWVfbWluaW1hbCgpCgpwcnNpY2VfZmVtYWxlX21ldGEgPC0gZnJlYWQoInByc2ljZV9uYWxsc19mZW1hbGVfZGF0YS9URVJSRV9QUlNpY2VfbmFsbHNfZmVtYWxlLnByc2ljZSIpCmdncGxvdChwcnNpY2VfZmVtYWxlX21ldGFbVGhyZXNob2xkIDw9IDAuNV0sIGFlcyhUaHJlc2hvbGQsIE51bV9TTlAsIGNvbG9yID0gLWxvZzEwKFApKSkgKwogIGdlb21fcG9pbnQoKSArCiAgc2NhbGVfeV9jb250aW51b3VzKGJyZWFrcyA9IGMoc2VxKDAsIDFlNSwgMi41ZTQpLCBzZXEoMmU1LCA2ZTUsIDFlNSkpKSArCiAgdGhlbWVfbWluaW1hbCgpCgpwcnNpY2VfbWV0YSA8LSBmcmVhZCgicHJzaWNlX2RhdGEvVEVSUkVfUFJTaWNlLnByc2ljZSIpCmdncGxvdChtYXBwaW5nID0gYWVzKFRocmVzaG9sZCwgUjIsIGNvbG9yID0gLWxvZzEwKFApKSkgKwogIGdlb21fcG9pbnQoZGF0YSA9IHByc2ljZV9mZW1hbGVfbWV0YSwgc2l6ZSA9IDEpICsKICBzY2FsZV9jb2xvcl9ncmFkaWVudChsb3cgPSAibGlnaHRwaW5rNCIsIGhpZ2ggPSAibGlnaHRwaW5rIikgKwogIGxhYnMoY29sb3IgPSBicXVvdGUoIkZlbWFsZSBsb2ciWyIxMCJdIH4gIihQKSIpKSArCiAgbmV3X3NjYWxlX2NvbG9yKCkgKwogIGdlb21fcG9pbnQoZGF0YSA9IHByc2ljZV9tYWxlX21ldGEsIHNpemUgPSAxLCBhZXMoY29sb3IgPSAtbG9nMTAoUCkpKSArCiAgc2NhbGVfY29sb3JfZ3JhZGllbnQobG93ID0gImxpZ2h0Ymx1ZTQiLCBoaWdoID0gImxpZ2h0Ymx1ZSIpICsKICBsYWJzKHkgPSBicXVvdGUoIlIiXjIpLCB4ID0gIkdXQVMgUC1WYWx1ZSBUaHJlc2hvbGQiLCBjb2xvciA9IGJxdW90ZSgiTWFsZSAtbG9nIlsiMTAiXSB+ICIoUCkiKSkgKwogIHRoZW1lX21pbmltYWwoKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gInRvcCIpCmdncGxvdChtYXBwaW5nID0gYWVzKFRocmVzaG9sZCwgUjIsIGNvbG9yID0gLWxvZzEwKFApKSkgKwogIGdlb21fcG9pbnQoZGF0YSA9IHByc2ljZV9mZW1hbGVfbWV0YSwgc2l6ZSA9IDEpICsKICBzY2FsZV9jb2xvcl9ncmFkaWVudChsb3cgPSAibGlnaHRwaW5rNCIsIGhpZ2ggPSAibGlnaHRwaW5rIixndWlkZSA9IGd1aWRlX2NvbG9yYmFyKG9yZGVyPTMpKSArCiAgbGFicyhjb2xvciA9IGJxdW90ZSgiRmVtYWxlIC1sb2ciWyIxMCJdIH4gIihQKSIpKSArCiAgbmV3X3NjYWxlX2NvbG9yKCkgKwogIGdlb21fcG9pbnQoZGF0YSA9IHByc2ljZV9tYWxlX21ldGEsIHNpemUgPSAxLCBhZXMoY29sb3IgPSAtbG9nMTAoUCkpKSArCiAgc2NhbGVfY29sb3JfZ3JhZGllbnQobG93ID0gImxpZ2h0Ymx1ZTQiLCBoaWdoID0gImxpZ2h0Ymx1ZSIsZ3VpZGUgPSBndWlkZV9jb2xvcmJhcihvcmRlcj0yKSkgKwogIGxhYnMoY29sb3IgPSBicXVvdGUoIk1hbGUgLWxvZyJbIjEwIl0gfiAiKFApIikpICsKICBuZXdfc2NhbGVfY29sb3IoKSArCiAgZ2VvbV9wb2ludChkYXRhID0gcHJzaWNlX21ldGEsIHNpemUgPSAxLCBhZXMoY29sb3IgPSAtbG9nMTAoUCkpKSArCiAgc2NhbGVfY29sb3JfZ3JhZGllbnQobG93ID0gImdyYXk0MCIsIGhpZ2ggPSAiZ3JheTgwIixndWlkZSA9IGd1aWRlX2NvbG9yYmFyKG9yZGVyPTEpKSArCiAgbGFicyh5ID0gYnF1b3RlKCJSIl4yKSwgeCA9ICJHV0FTIFAtVmFsdWUgVGhyZXNob2xkIiwgY29sb3IgPSBicXVvdGUoIkNyb3NzLXNleCAtbG9nIlsiMTAiXSB+ICIoUCkiKSkgKwogIHRoZW1lX21pbmltYWwoKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gInRvcCIsbGVnZW5kLnRpdGxlID0gZWxlbWVudF90ZXh0KHNpemU9NykpCmBgYAoKIyBTdGVwIDMgcnVuIGxpbmVhciBtb2RlbCBhdCBkaWZmZXJlbnQgdGhyZXNob2xkcyBmb3IgU05QIGluY2x1c2lvbgpUaGlzIGlzIHVwZGF0ZWQgdG8gcmVjZW50bHkgcHJvY2Vzc2VkIEROQW0gRGF0YSAyMCBvdXRsaWVycyBpbiBETkFtIHJlbW92ZWQ6CmBgYHtyfQpsb2FkKCIvaG9tZTEvTkVVUk8vU0hBUkVfREVDSVBIRVIvcHJvY2Vzc2VkX0ROQW1fZGF0YS8yMDIyL1RFUlJFX3Byb2Nlc3NlZF8yMDIyLzEtVEVSUkVfUkdfZmlsdGVyZWQuUkRhdGEiKSAjUERfUkdfZmlsdGVyZWQKIyBBc3NpZ24gZ2Vub3R5cGluZyBJRCB0byBkYXRhCm9yaWdpbmFsX2NvdmFycyA8LSBmcmVhZCgiL2hvbWUxL05FVVJPL1NIQVJFX0RFQ0lQSEVSL3RlcnJlX21ldGFfbWFzdGVyLmNzdiIpWywgLihwYXRpZW50LCBJSUQgPSBnc3ViKCJfUEFFLioiLCAiIiwgSUlEKSldCmJldGFzX2NvbWJhdCA8LSBtaW5maTo6Z2V0QmV0YShQRF9SR19maWx0ZXJlZCkKY29sbmFtZXMoYmV0YXNfY29tYmF0KSA8LSBvcmlnaW5hbF9jb3ZhcnMkSUlEW21hdGNoKGNvbG5hbWVzKGJldGFzX2NvbWJhdCksIG9yaWdpbmFsX2NvdmFycyRwYXRpZW50KV0KYmV0YXNfY29tYmF0IDwtIGJldGFzX2NvbWJhdFssIGNvbG5hbWVzKGJldGFzX2NvbWJhdCkgJWluJSBjb3ZhcmlhdGUkSUlEXQpgYGAKTGV0J3MgY2hlY2sgaG93IHRoZSBkYXRhIGxvb2tzIGZvciB0aGUgZmlyc3QgNSBzdWJqZWN0czoKYGBge3J9CmdncGxvdChiZXRhc19jb21iYXRbLCAxOjVdICU+JSBhcy5kYXRhLnRhYmxlKGtlZXAucm93bmFtZXMgPSBUKSAlPiUgbWVsdChpZC52YXJzID0gInJuIiwgdmFsdWUubmFtZSA9ICJiZXRhcyIsIHZhcmlhYmxlLm5hbWUgPSAic3ViamVjdCIpLCBhZXMoYmV0YXMsIGNvbG9yID0gc3ViamVjdCkpKwogIGdlb21fZGVuc2l0eSgpCmBgYAoKIyMjIE1hdGNoIEROQSwgUFJTLCBhbmQgbWV0YWRhdGEgQFRPRE8gZml4IHRoaXMgYWZ0ZXIgdXBkYXRlIG9mIFBSUyBldGMuCmBgYHtyfQojIHByc2ljZV9iZXN0PTAuMDAwMzAwMDUKcHJzaWNlX2FsbCA8LSBmcmVhZCgicHJzaWNlX2RhdGEvVEVSUkVfUFJTaWNlLmFsbF9zY29yZSIpW21hdGNoKGNvbG5hbWVzKGJldGFzX2NvbWJhdCksIElJRCxub21hdGNoPTApLCAuKEZJRCwgSUlELCBgUHRfNWUtMDhgLCBgUHRfNS4wMDVlLTA1YCwgYFB0XzAuMDAwMTAwMDVgLCBgUHRfMC4wMDEwMDAwNWAsIGBQdF8wLjAwMDMwMDA1YCxgUHRfMC4wMTAxNTAxYCwgYFB0XzAuMWAsIGBQdF8wLjJgLCBgUHRfMC4zYCwgYFB0XzAuNGAsIGBQdF8wLjVgLCBgUHRfMWApXQpjb3ZhcmlhdGUgPC0gY292YXJpYXRlW21hdGNoKGNvbG5hbWVzKGJldGFzX2NvbWJhdCksIElJRCldCmFsbChjb3ZhcmlhdGUkSUlEID09IGNvbG5hbWVzKGJldGFzX2NvbWJhdCkpCmFsbChjb3ZhcmlhdGUkSUlEID09IHByc2ljZV9hbGwkSUlEKQoKY292YXJpYXRlX21hbGUgPC0gY292YXJpYXRlW3NleCA9PSAxXSAlPiUgc2VsZWN0KC1zZXgpCmJldGFzX21hbGUgPC0gYmV0YXNfY29tYmF0WywgY292YXJpYXRlX21hbGUkSUlEXQojIGJlc3QgPSAwLjAzMzUwMDEKcHJzaWNlX21hbGVfYWxsIDwtZnJlYWQoInByc2ljZV9uYWxsc19tYWxlX2RhdGEvVEVSUkVfUFJTaWNlX25hbGxzX21hbGUuYWxsX3Njb3JlIilbbWF0Y2goY29sbmFtZXMoYmV0YXNfbWFsZSksIElJRCxub21hdGNoPTApLCAuKEZJRCwgSUlELCBgUHRfNWUtMDhgLCBgUHRfNS4wMDVlLTA1YCwgYFB0XzAuMDAwMTAwMDVgLCBgUHRfMC4wMDEwMDAwNWAsIGBQdF8wLjAxMDE1MDFgLCBgUHRfMC4wMzM1MDAxYCwgYFB0XzAuMWAsIGBQdF8wLjJgLCBgUHRfMC4zYCwgYFB0XzAuNGAsIGBQdF8wLjVgLCBgUHRfMWApXQoKY292YXJpYXRlX21hbGUgPC0gY292YXJpYXRlX21hbGVbbWF0Y2goY29sbmFtZXMoYmV0YXNfbWFsZSksIElJRCxub21hdGNoPTApXQphbGwoY292YXJpYXRlX21hbGUkSUlEID09IGNvbG5hbWVzKGJldGFzX21hbGUpKQphbGwoY292YXJpYXRlX21hbGUkSUlEID09IHByc2ljZV9tYWxlX2FsbCRJSUQpCgpjb3ZhcmlhdGVfZmVtYWxlIDwtIGNvdmFyaWF0ZVtzZXggPT0gMF0gJT4lIHNlbGVjdCgtc2V4KSAjJT4lIGZpbHRlcihJSUQgJWluJSBwcnNpY2VfZmVtYWxlX2Jlc3QkSUlEKQpiZXRhc19mZW1hbGUgPC0gYmV0YXNfY29tYmF0WywgY292YXJpYXRlX2ZlbWFsZSRJSURdCiNiZXN0ID0gMC4wMDAzNTAwNQoKcHJzaWNlX2ZlbWFsZV9hbGwgPC0gZnJlYWQoInByc2ljZV9uYWxsc19mZW1hbGVfZGF0YS9URVJSRV9QUlNpY2VfbmFsbHNfZmVtYWxlLmFsbF9zY29yZSIpW21hdGNoKGNvbG5hbWVzKGJldGFzX2ZlbWFsZSksIElJRCxub21hdGNoPTApLCAuKEZJRCwgSUlELCBgUHRfNWUtMDhgLCBgUHRfNS4wMDVlLTA1YCwgYFB0XzAuMDAwMTAwMDVgLCBgUHRfMC4wMDAzNTAwNWAsIGBQdF8wLjAwMTAwMDA1YCwgYFB0XzAuMDEwMTUwMWAsIGBQdF8wLjFgLCBgUHRfMC4yYCwgYFB0XzAuM2AsIGBQdF8wLjRgLCBgUHRfMC41YCwgYFB0XzFgKV0KY292YXJpYXRlX2ZlbWFsZSA8LSBjb3ZhcmlhdGVfZmVtYWxlW21hdGNoKGNvbG5hbWVzKGJldGFzX2ZlbWFsZSksIElJRCxub21hdGNoPTApXQphbGwoY292YXJpYXRlX2ZlbWFsZSRJSUQgPT0gY29sbmFtZXMoYmV0YXNfZmVtYWxlKSkKYWxsKGNvdmFyaWF0ZV9mZW1hbGUkSUlEID09IHByc2ljZV9mZW1hbGVfYWxsJElJRCkKCgpgYGAKIyMjIFJ1biBsaW1tYQoKYGBge3J9Cm12YWx1ZXMgPC0gbHVtaTo6YmV0YTJtKGJldGFzX2NvbWJhdCkKcHJzX21hdCA8LSBwcnNpY2VfYWxsWywgLWMoMSwgMildCmNvdl9tYXQgPC0gY292YXJpYXRlWywgLWMoMSwgMildCgptdmFsdWVzX21hbGUgPC0gbHVtaTo6YmV0YTJtKGJldGFzX21hbGUpCnByc19tYXRfbWFsZSA8LSBwcnNpY2VfbWFsZV9hbGxbLCAtYygxLCAyKV0KY292X21hdF9tYWxlIDwtIGNvdmFyaWF0ZV9tYWxlWywgLWMoMSwgMildCgptdmFsdWVzX2ZlbWFsZSA8LSBsdW1pOjpiZXRhMm0oYmV0YXNfZmVtYWxlKQpwcnNfbWF0X2ZlbWFsZSA8LSBwcnNpY2VfZmVtYWxlX2FsbFssIC1jKDEsIDIpXQpjb3ZfbWF0X2ZlbWFsZSA8LSBjb3ZhcmlhdGVfZmVtYWxlWywgLWMoMSwgMildCmBgYAoKYGBge3J9CnJlZ2lzdGVyRG9QYXJhbGxlbChuY29sKHByc19tYXQpIC8gNCkKaGl0cyA8LSBmb3JlYWNoKHByc190aHJlc2ggPSBjb2xuYW1lcyhwcnNfbWF0KSkgJWRvcGFyJSB7CiAgZGVzaWduX3BycyA8LSBtb2RlbC5tYXRyaXgofi4sIGRhdGEgPSBjYmluZChwcnNfbWF0WywgLi5wcnNfdGhyZXNoXSwgY292X21hdCkpCiAgcHJzX2ZpdCA8LSBsbUZpdChtdmFsdWVzLCBkZXNpZ25fcHJzKQogIHByc19maXQgPC0gZUJheWVzKHByc19maXQpCiAgdG9wVGFibGUocHJzX2ZpdCwgY29lZiA9IDIsIGFkanVzdC5tZXRob2QgPSAiYm9uZiIsIHAudmFsdWUgPSAwLjA1LCBudW1iZXIgPSBJbmYsIGdlbmVsaXN0ID0gcm93bmFtZXMobXZhbHVlcykpCn0KbmFtZXMoaGl0cykgPC0gY29sbmFtZXMocHJzX21hdCkKaGl0c19ieV90aHJlc2hfYm9uZiA8LSByYmluZGxpc3QoaGl0cywgaWRjb2wgPSAidGhyZXNob2xkIiwgZmlsbCA9IFRSVUUpCgpyZWdpc3RlckRvUGFyYWxsZWwobmNvbChwcnNfbWF0X21hbGUpIC8gNCkKaGl0c19tYWxlIDwtIGZvcmVhY2gocHJzX3RocmVzaCA9IGNvbG5hbWVzKHByc19tYXRfbWFsZSkpICVkb3BhciUgewogIGRlc2lnbl9wcnNfbWFsZSA8LSBtb2RlbC5tYXRyaXgofi4sIGRhdGEgPSBjYmluZChwcnNfbWF0X21hbGVbLCAuLnByc190aHJlc2hdLCBjb3ZfbWF0X21hbGUpKQogIHByc19maXRfbWFsZSA8LSBsbUZpdChtdmFsdWVzX21hbGUsIGRlc2lnbl9wcnNfbWFsZSkKICBwcnNfZml0X21hbGUgPC0gZUJheWVzKHByc19maXRfbWFsZSkKICB0b3BUYWJsZShwcnNfZml0X21hbGUsIGNvZWYgPSAyLCBhZGp1c3QubWV0aG9kID0gImJvbmYiLCBwLnZhbHVlID0gMC4wNSwgbnVtYmVyID0gSW5mLCBnZW5lbGlzdCA9IHJvd25hbWVzKG12YWx1ZXNfbWFsZSkpCn0KbmFtZXMoaGl0c19tYWxlKSA8LSBjb2xuYW1lcyhwcnNfbWF0X21hbGUpCmhpdHNfYnlfdGhyZXNoX2JvbmZfbWFsZSA8LSByYmluZGxpc3QoaGl0c19tYWxlLCBpZGNvbCA9ICJ0aHJlc2hvbGQiLCBmaWxsID0gVFJVRSkKCnJlZ2lzdGVyRG9QYXJhbGxlbChuY29sKHByc19tYXRfZmVtYWxlKSAvIDQpCmhpdHNfZmVtYWxlIDwtIGZvcmVhY2gocHJzX3RocmVzaCA9IGNvbG5hbWVzKHByc19tYXRfZmVtYWxlKSkgJWRvcGFyJSB7CiAgZGVzaWduX3Byc19mZW1hbGUgPC0gbW9kZWwubWF0cml4KH4uLCBkYXRhID0gY2JpbmQocHJzX21hdF9mZW1hbGVbLCAuLnByc190aHJlc2hdLCBjb3ZfbWF0X2ZlbWFsZSkpCiAgcHJzX2ZpdF9mZW1hbGUgPC0gbG1GaXQobXZhbHVlc19mZW1hbGUsIGRlc2lnbl9wcnNfZmVtYWxlKQogIHByc19maXRfZmVtYWxlIDwtIGVCYXllcyhwcnNfZml0X2ZlbWFsZSkKICB0b3BUYWJsZShwcnNfZml0X2ZlbWFsZSwgY29lZiA9IDIsIGFkanVzdC5tZXRob2QgPSAiYm9uZiIsIHAudmFsdWUgPSAwLjA1LCBudW1iZXIgPSBJbmYsIGdlbmVsaXN0ID0gcm93bmFtZXMobXZhbHVlc19mZW1hbGUpKQp9Cm5hbWVzKGhpdHNfZmVtYWxlKSA8LSBjb2xuYW1lcyhwcnNfbWF0X2ZlbWFsZSkKaGl0c19ieV90aHJlc2hfYm9uZl9mZW1hbGUgPC0gcmJpbmRsaXN0KGhpdHNfZmVtYWxlLCBpZGNvbCA9ICJ0aHJlc2hvbGQiLCBmaWxsID0gVFJVRSkKCmBgYAoKIyMjIFBsb3R0aW5nIEVXQVMgdnMgVGhyZXNob2xkIEV4cGVyaW1lbnQgYnkgU2V4CmBgYHtyfQp0b19wbG90IDwtIHJiaW5kKAogIGhpdHNfYnlfdGhyZXNoX2JvbmZbLCAuKGhpdHMgPSAuTiwgU2V4ID0gIkNyb3NzLXNleCIpLCBieSA9IHRocmVzaG9sZF0gJT4lCiAgICBtdXRhdGUodGhyZXNob2xkID0gcmVjb2RlX2ZhY3Rvcih0aHJlc2hvbGQsIGBQdF8wLjAyMTkwMDFgID0gIjAuMDIxOSIsIGBQdF81ZS0wOGAgPSAiNWUtOCIsIGBQdF81LjAwNWUtMDVgID0gIjVlLTUiLCBgUHRfMC4wMDAxMDAwNWAgPSAiMWUtNCIsYFB0XzAuMDAwMzAwMDVgPSIzZS00IixgUHRfMC4wMDEwMDAwNWAgPSAiMWUtMyIsIGBQdF8wLjAxMDE1MDFgID0gIjFlLTIiLCBgUHRfMC4xYCA9ICIwLjEiLCBgUHRfMC4yYCA9ICIwLjIiLCBgUHRfMC4zYCA9ICIwLjMiLCBgUHRfMC40YCA9ICIwLjQiLCBgUHRfMC41YCA9ICIwLjUiLCBgUHRfMWAgPSAiMS4wIikpLAogIGhpdHNfYnlfdGhyZXNoX2JvbmZfbWFsZVssIC4oaGl0cyA9IC5OLCBTZXggPSAiTWFsZSIpLCBieSA9IHRocmVzaG9sZF0gJT4lCiAgICBtdXRhdGUodGhyZXNob2xkID0gcmVjb2RlX2ZhY3Rvcih0aHJlc2hvbGQsIGBQdF8wLjAyMTkwMDFgID0gIjAuMDIxOSIsIGBQdF81ZS0wOGAgPSAiNWUtOCIsIGBQdF81LjAwNWUtMDVgID0gIjVlLTUiLCBgUHRfMC4wMDAxMDAwNWAgPSAiMWUtNCIsIGBQdF8wLjAwMTAwMDA1YCA9ICIxZS0zIiwgYFB0XzAuMDEwMTUwMWAgPSAiMWUtMiIsIGBQdF8wLjAzMzUwMDFgPSIzLjRlLTIiLCBgUHRfMC4xYCA9ICIwLjEiLCBgUHRfMC4yYCA9ICIwLjIiLCBgUHRfMC4zYCA9ICIwLjMiLCBgUHRfMC40YCA9ICIwLjQiLCBgUHRfMC41YCA9ICIwLjUiLCBgUHRfMWAgPSAiMS4wIikpLAogIGhpdHNfYnlfdGhyZXNoX2JvbmZfZmVtYWxlWywgLihoaXRzID0gLk4sIFNleCA9ICJGZW1hbGUiKSwgYnkgPSB0aHJlc2hvbGRdICU+JQogICAgbXV0YXRlKHRocmVzaG9sZCA9IHJlY29kZV9mYWN0b3IodGhyZXNob2xkLCBgUHRfMC4wMjE5MDAxYCA9ICIwLjAyMTkiLCBgUHRfNWUtMDhgID0gIjVlLTgiLCBgUHRfNS4wMDVlLTA1YCA9ICI1ZS01IiwgYFB0XzAuMDAwMTAwMDVgID0gIjFlLTQiLGBQdF8wLjAwMDM1MDA1YD0iMy41ZS00IiwgYFB0XzAuMDAxMDAwMDVgID0gIjFlLTMiLCBgUHRfMC4wMTAxNTAxYCA9ICIxZS0yIiwgYFB0XzAuMWAgPSAiMC4xIiwgYFB0XzAuMmAgPSAiMC4yIiwgYFB0XzAuM2AgPSAiMC4zIiwgYFB0XzAuNGAgPSAiMC40IiwgYFB0XzAuNWAgPSAiMC41IiwgYFB0XzFgID0gIjEuMCIpKQopICU+JSBtdXRhdGUoU2V4ID0gZmFjdG9yKFNleCwgbGV2ZWxzID0gYygiQ3Jvc3Mtc2V4IiwgIk1hbGUiLCAiRmVtYWxlIikpKQpwbG90X3BvcyA8LXBvc2l0aW9uX2RvZGdlMih3aWR0aCA9IDAuOSwgcHJlc2VydmUgPSAic2luZ2xlIikKZ2dwbG90KHRvX3Bsb3QsIGFlcyh0aHJlc2hvbGQsIGhpdHMsIGZpbGwgPSBTZXgsIGxhYmVsID0gaGl0cykpICsKICBnZW9tX3RleHQocG9zaXRpb24gPSBwbG90X3Bvcywgdmp1c3QgPSAtMC4yNSkgKwogIGdlb21fY29sKHBvc2l0aW9uID0gcGxvdF9wb3MpICsKICBsYWJzKHggPSAiR1dBUyBQIFZhbHVlIFRocmVzaG9sZCIsIHkgPSAiRVdBUyBIaXRzIikgKwogIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcyA9IGMoImdyZXk4MCIsICJsaWdodGJsdWUiLCAibGlnaHRwaW5rIikpICsKICB0aGVtZV9taW5pbWFsKCkKaGl0c19ieV90aHJlc2hfYm9uZlssIC4oaGl0cyA9IC5OLCBTZXggPSAiQ3Jvc3Mtc2V4IiksIGJ5ID0gdGhyZXNob2xkXQpoaXRzX2J5X3RocmVzaF9ib25mX21hbGVbLCAuKGhpdHMgPSAuTiksIGJ5ID0gdGhyZXNob2xkXQpoaXRzX2J5X3RocmVzaF9ib25mX2ZlbWFsZVssIC4oaGl0cyA9IC5OKSwgYnkgPSB0aHJlc2hvbGRdCmBgYApgYGB7cn0KZGlzcGxheV92ZW5uIDwtIGZ1bmN0aW9uKHgsIC4uLikgewogIGxpYnJhcnkoVmVubkRpYWdyYW0pCiAgZ3JpZC5uZXdwYWdlKCkKICB2ZW5uX29iamVjdCA8LSB2ZW5uLmRpYWdyYW0oeCwgZmlsZW5hbWUgPSBOVUxMLCAuLi4pCiAgZ3JpZC5kcmF3KHZlbm5fb2JqZWN0KQp9CgpkaXNwbGF5X3Zlbm4obGlzdChgQ3Jvc3Mtc2V4YCA9IGhpdHNfYnlfdGhyZXNoX2JvbmZbdGhyZXNob2xkID09ICJQdF81ZS0wOCJdJElELCBNYWxlID0gaGl0c19ieV90aHJlc2hfYm9uZl9tYWxlW3RocmVzaG9sZCA9PSAiUHRfNWUtMDgiXSRJRCwgRmVtYWxlID0gaGl0c19ieV90aHJlc2hfYm9uZl9mZW1hbGVbdGhyZXNob2xkID09ICJQdF81ZS0wOCJdJElEKSwgZmlsbCA9IGMoImdyYXk4MCIsICJsaWdodGJsdWUiLCAibGlnaHRwaW5rIikpCm1hbGVfaWRzIDwtIGhpdHNfYnlfdGhyZXNoX2JvbmZfbWFsZVt0aHJlc2hvbGQgPT0gIlB0XzVlLTA4Il0kSUQKY3Jvc3NfaWRzIDwtIGhpdHNfYnlfdGhyZXNoX2JvbmZbdGhyZXNob2xkID09ICJQdF81ZS0wOCJdJElECm1hbGVfaWRzWyFtYWxlX2lkcyAlaW4lIGNyb3NzX2lkc10KYGBgCgpgYGB7cn0KZ2V0X2Z1bGxfZml0IDwtIGZ1bmN0aW9uKHByc19tYXQsY292X21hdCxtdmFsdWVzKXsKICB0b3BfZGVzaWduX3BycyA8LSBtb2RlbC5tYXRyaXgofi4sIGRhdGEgPSBjYmluZChwcnNfbWF0WywgYFB0XzVlLTA4YF0sIGNvdl9tYXQpKQogIHRvcF9wcnNfZml0IDwtIGxtRml0KG12YWx1ZXMsIHRvcF9kZXNpZ25fcHJzKQogIHRvcF9wcnNfZml0IDwtIGVCYXllcyh0b3BfcHJzX2ZpdCkKICB0b3BfcHJzX2hpdHMgPC0gdG9wVGFibGUodG9wX3Byc19maXQsIGNvZWYgPSAyLCBhZGp1c3QubWV0aG9kID0gImJvbmYiLCBudW1iZXIgPSBJbmYsIGdlbmVsaXN0ID0gcm93bmFtZXMobXZhbHVlcykpCn0KdG9wX3Byc19oaXRzIDwtIGdldF9mdWxsX2ZpdChwcnNfbWF0LGNvdl9tYXQsbXZhbHVlcykKdG9wX21hbGVfcHJzX2hpdHMgPC0gZ2V0X2Z1bGxfZml0KHByc19tYXRfbWFsZSwgY292X21hdF9tYWxlLCBtdmFsdWVzX21hbGUpCnRvcF9mZW1hbGVfcHJzX2hpdHMgPC0gZ2V0X2Z1bGxfZml0KHByc19tYXRfZmVtYWxlLCBjb3ZfbWF0X2ZlbWFsZSwgbXZhbHVlc19mZW1hbGUpCnNhdmUobGlzdD1jKCJ0b3BfcHJzX2hpdHMiLCJ0b3BfbWFsZV9wcnNfaGl0cyIsInRvcF9mZW1hbGVfcHJzX2hpdHMiKSxmaWxlPSJwcnNfbmFsbHNfY3Jvc3Nfd19zZXhfc3RyYXRpZmllZC5SRGF0YSIpCmBgYAoKYGBge3J9CmxvYWQoInByc19uYWxsc19jcm9zc193X3NleF9zdHJhdGlmaWVkLlJEYXRhIikKYGBgCmBgYHtyfQptYW5pZmVzdCA8LSBJbGx1bWluYUh1bWFuTWV0aHlsYXRpb25FUElDYW5uby5pbG0xMGI0LmhnMTk6Ok90aGVyICU+JQogIGFzLmRhdGEuZnJhbWUoKSAlPiUKICByb3duYW1lc190b19jb2x1bW4odmFyID0gIm5hbWUiKQpwcnNfYW5ub3QgPC0gZGF0YS50YWJsZSh0b3BfcHJzX2hpdHMpW21hbmlmZXN0LCBnZW5lIDo9IGdzdWIoIjsuKiIsICIiLCBVQ1NDX1JlZkdlbmVfTmFtZSksIG9uID0gYyhJRCA9ICJuYW1lIildCnByc19hbm5vdF9tYWxlIDwtIGRhdGEudGFibGUodG9wX21hbGVfcHJzX2hpdHMpW21hbmlmZXN0LCBnZW5lIDo9IGdzdWIoIjsuKiIsICIiLCBVQ1NDX1JlZkdlbmVfTmFtZSksIG9uID0gYyhJRCA9ICJuYW1lIildCnByc19hbm5vdF9mZW1hbGU8LSBkYXRhLnRhYmxlKHRvcF9mZW1hbGVfcHJzX2hpdHMpW21hbmlmZXN0LCBnZW5lIDo9IGdzdWIoIjsuKiIsICIiLCBVQ1NDX1JlZkdlbmVfTmFtZSksIG9uID0gYyhJRCA9ICJuYW1lIildCnBsb3RfcHJzX2hpdHMgPC0gZnVuY3Rpb24ocHJzX2Fubm90LGxhYmVsX2NvbG9yKXsKICBnZ3Bsb3QocHJzX2Fubm90LCBhZXMobG9nRkMsIC1sb2cxMChQLlZhbHVlKSkpICsKICAgIGdlb21fcG9pbnQoKSArCiAgICBnZW9tX3BvaW50KAogICAgICBkYXRhID0gc3Vic2V0KHByc19hbm5vdCwgYWRqLlAuVmFsIDwgMC4wNSAmIGFicyhsb2dGQykgPiAwLjAzKSwKICAgICAgY29sb3IgPSBsYWJlbF9jb2xvciwKICAgICAgbWFwcGluZyA9IGFlcyhsb2dGQywgLWxvZzEwKFAuVmFsdWUpKQogICAgKSArCiAgICBnZW9tX2hsaW5lKAogICAgICBsaW5ldHlwZSA9ICJkYXNoZWQiLAogICAgICB5aW50ZXJjZXB0ID0gbWluKC1sb2cxMChwcnNfYW5ub3QkUC5WYWx1ZVtwcnNfYW5ub3QkYWRqLlAuVmFsIDwgMC4wNV0pKQogICAgKSArCiAgICBnZW9tX3ZsaW5lKGxpbmV0eXBlID0gImRhc2hlZCIsIHhpbnRlcmNlcHQgPSAwLjAzKSArCiAgICBnZW9tX3ZsaW5lKGxpbmV0eXBlID0gImRhc2hlZCIsIHhpbnRlcmNlcHQgPSAtMC4wMykgKwogICAgZ2VvbV90ZXh0X3JlcGVsKAogICAgICBkYXRhID0gcHJzX2Fubm90ICU+JSBmaWx0ZXIoYWJzKGxvZ0ZDKSA+IDAuMDMgJiBhZGouUC5WYWwgPCAwLjA1KSwKICAgICAgY29sb3IgPSAiZG9kZ2VyYmx1ZSIsCiAgICAgIG1hcHBpbmcgPSBhZXMobG9nRkMsIC1sb2cxMChQLlZhbHVlKSwgbGFiZWwgPSBpZmVsc2UoZ2VuZSAhPSAiIiwgZ2VuZSwgSUQpKSwKICAgICAgc2l6ZSA9IDMsCiAgICAgIG1heC5vdmVybGFwcyA9IDIwCiAgICApICsKICAgIGxhYnMoeSA9IGJxdW90ZSgibG9nIlsxMF0gfiAiKFApIiksIHggPSBxdW90ZShEZWx0YSB+ICJNIiB+IE1ldGh5bGF0aW9uKSkgKwogICAgdGhlbWVfbWluaW1hbCgpCn0KCnBsb3RfcHJzX2hpdHMocHJzX2Fubm90LCJncmF5NDAiKQpwbG90X3Byc19oaXRzKHByc19hbm5vdF9tYWxlLCJsaWdodGJsdWUiKQpwbG90X3Byc19oaXRzKHByc19hbm5vdF9mZW1hbGUsInBpbmsiKQpgYGAKCgpgYGB7cn0KcHJzX3Bsb3RfZGF0YSA8LSBkYXRhLnRhYmxlKGRuYW0gPSBsdW1pOjptMmJldGEobXZhbHVlc1siY2cxMjYwOTc4NSIsXSksIHBycyA9IHByc19tYXQkYFB0XzVlLTA4YCxQRD1pZmVsc2UoUERbbWF0Y2goY29sbmFtZXMobXZhbHVlcyksSUlEKSxdJFBEID09IDEsIkNhc2UiLCJDb250cm9sIiksU2V4PWZhY3RvcihpZmVsc2UoY292X21hdCRzZXggPT0xLCJNYWxlIiwiRmVtYWxlIiksbGV2ZWxzPWMoIk1hbGUiLCJGZW1hbGUiKSkpCmdncGxvdChwcnNfcGxvdF9kYXRhLGFlcyhwcnMsZG5hbSxjb2xvcj1TZXgsc2hhcGUgPSBQRCxncm91cD1TZXgpKSArCiAgZ2VvbV9wb2ludChzaXplPTMpKwogIGdlb21fc21vb3RoKG1ldGhvZD0ibG0iLHNlPUYpKwogIGxhYnMoeT1icXVvdGUoIk1ldGh5bGF0aW9uIn5iZXRhKSx4PSJQYXJraW5zb24ncyBHUlMiKSsKICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzPWMoIk1hbGUiPSJsaWdodGJsdWUiLCJGZW1hbGUiPSJsaWdodHBpbmsiKSkgKyAKICB0aGVtZV9taW5pbWFsKGJhc2Vfc2l6ZT0yMCkKYGBgCgojIyBnbyBlbnJpY2htZW50PwpgYGB7cixldmFsPUZBTFNFfQpsaWJyYXJ5KGdwcm9maWxlcjIpCmNyb3NzX3NleCA8LSB1bmlxdWUoZ3N1YigiOy4qIiwiIixtYW5pZmVzdFttYW5pZmVzdCRuYW1lICVpbiUgaGl0c19ieV90aHJlc2hfYm9uZlt0aHJlc2hvbGQgPT0gIlB0XzVlLTA4Il0kSUQsXSRVQ1NDX1JlZkdlbmVfTmFtZSkpCmJhY2tncm91bmQgPC0gdW5pcXVlKGdzdWIoIjsuKiIsIiIsbWFuaWZlc3QkVUNTQ19SZWZHZW5lX05hbWUpKQpnb3N0X3JlcyA8LSBnb3N0KHF1ZXJ5PWNyb3NzX3NleCxjdXN0b21fYmcgPSBiYWNrZ3JvdW5kKQpnb3N0cGxvdChnb3N0X3JlcykKYGBgCmBgYHtyfQpjcGdfcG9zIDwtIElsbHVtaW5hSHVtYW5NZXRoeWxhdGlvbkVQSUNhbm5vLmlsbTEwYjQuaGcxOTo6TG9jYXRpb25zICU+JQogIGFzLmRhdGEuZnJhbWUoKSAlPiUKICByb3duYW1lc190b19jb2x1bW4odmFyID0gIm5hbWUiKQptYWxlX2Fubm90X2NwZyA8LSAocHJzX2Fubm90X21hbGUgJT4lIGZpbHRlcihhZGouUC5WYWwgPCAwLjI1KSAlPiUgbGVmdF9qb2luKGNwZ19wb3MsYnk9YygiSUQiPSJuYW1lIikpICU+JSBncm91cF9ieShjaHIpICU+JSB0b3BfbigtMTAsUC5WYWx1ZSkpJElECmZlbWFsZV9hbm5vdF9jcGcgPC0gKHByc19hbm5vdF9mZW1hbGUgJT4lIGZpbHRlcihhZGouUC5WYWwgPCAwLjI1KSAlPiUgbGVmdF9qb2luKGNwZ19wb3MsYnk9YygiSUQiPSJuYW1lIikpICU+JSBncm91cF9ieShjaHIpICU+JSB0b3BfbigtMTAsUC5WYWx1ZSkpJElECmxpYnJhcnkoaHVkc29uKQpvcHRpb25zKGdncmVwZWwubWF4Lm92ZXJsYXBzID0gSW5mKQpnbWlycm9yKAogIHByc19hbm5vdF9tYWxlWywuKFNOUD1JRCxDSFI9Z3N1YigiY2hyIiwiIixjcGdfcG9zW21hdGNoKElELGNwZ19wb3MkbmFtZSksXSRjaHIpLFBPUz1jcGdfcG9zW21hdGNoKElELGNwZ19wb3MkbmFtZSksXSRwb3MscHZhbHVlPVAuVmFsdWUpXSwKICBwcnNfYW5ub3RfZmVtYWxlWywuKFNOUD1JRCxDSFI9Z3N1YigiY2hyIiwiIixjcGdfcG9zW21hdGNoKElELGNwZ19wb3MkbmFtZSksXSRjaHIpLFBPUz1jcGdfcG9zW21hdGNoKElELGNwZ19wb3MkbmFtZSksXSRwb3MscHZhbHVlPVAuVmFsdWUpXSwKICBhbm5vdGF0ZV9zbnAgPSBjKG1hbGVfYW5ub3RfY3BnLGZlbWFsZV9hbm5vdF9jcGcpLAogIHRsaW5lID0gIG1heChwcnNfYW5ub3RfbWFsZVthZGouUC5WYWw8MC4yNV0kUC5WYWx1ZSksCiAgYmxpbmUgPSAgbWF4KHByc19hbm5vdF9mZW1hbGVbYWRqLlAuVmFsPDAuMjVdJFAuVmFsdWUpLAogIGhpZ2hsaWdodF9wID0gYyhtYXgocHJzX2Fubm90X21hbGVbYWRqLlAuVmFsPDAuMjVdJFAuVmFsdWUpLG1heChwcnNfYW5ub3RfZmVtYWxlW2Fkai5QLlZhbDwwLjI1XSRQLlZhbHVlKSksCiAgdG9wdGl0bGU9Ik1hbGUiLAogIGJvdHRvbXRpdGxlID0iRmVtYWxlIiwKICBoaWdobGlnaHRlcj0iZ3JlZW4iLAogIGJhY2tncm91bmQ9IndoaXRlIgopCgpgYGAKCiMgTWFuaGF0dGFuIHBsb3QgdnMgR1dBUyBtYW5oYXR0YW4gcGxvdApgYGB7cixldmFsPUZBTFNFfQpsaWJyYXJ5KHFxbWFuKQoKY29weV9hbm5vdCA8LSBwcnNfYW5ub3RbY3BnX3Bvcywgb24gPSBjKElEID0gIm5hbWUiKV0gJT4lCiAgbXV0YXRlKGNociA9IGFzLm51bWVyaWMocmVjb2RlKGdzdWIoImNociIsICIiLCBjaHIpLFg9IjIzIixZPSIyNCIpKSkKdG9fcGxvdCA8LSBjb3B5X2Fubm90WywgLihTTlAgPSBnZW5lLCBDSFIgPSBjaHIsIEJQID0gcG9zLCBQID0gUC5WYWx1ZSwgRkRSID0gYWRqLlAuVmFsKV1bIWlzLm5hKFApXQptYW5oYXR0YW4odG9fcGxvdCwKICBhbm5vdGF0ZVB2YWwgPSBtYXgodG9fcGxvdFtGRFIgPCAwLjA1XSRQKSwKICBhbm5vdGF0ZVRvcCA9IFRSVUUsCiAgZ2Vub21ld2lkZWxpbmUgPSBtaW4oLWxvZzEwKHRvX3Bsb3RbRkRSIDwgMC4wNV0kUCkpLAogIHN1Z2dlc3RpdmVsaW5lID0gRkFMU0UKKQptYW5oYXR0YW4odG9fcGxvdFtDSFIgPT0gMTddLAogIGFubm90YXRlUHZhbCA9IG1heCh0b19wbG90W0ZEUiA8IDAuMDVdJFApLAogIGFubm90YXRlVG9wID0gRkFMU0UsCiAgZ2Vub21ld2lkZWxpbmUgPSBtaW4oLWxvZzEwKHRvX3Bsb3RbRkRSIDwgMC4wNV0kUCkpLAogIHN1Z2dlc3RpdmVsaW5lID0gRkFMU0UKKQpxcSh0b19wbG90JFApCmBgYAoKCiMjIERNUnMKCmBgYHtyfQpsaWJyYXJ5KERNUmNhdGUpClM0X3RvX2RhdGFmcmFtZSA8LSBmdW5jdGlvbihzNG9iaikgewogIG5tcyA8LSBzbG90TmFtZXMoczRvYmopCiAgbHN0IDwtIGxhcHBseShubXMsIGZ1bmN0aW9uKG5tKSBzbG90KHM0b2JqLCBubSkpCiAgYXMuZGF0YS5mcmFtZShzZXROYW1lcyhsc3QsIG5tcykpCn0KcnVuX2RtcmNhdGUgPC0gZnVuY3Rpb24ocHJzX21hdCxjb3ZfbWF0LG12YWx1ZXMpewogIGRlc2lnbl9wcnMgPC0gbW9kZWwubWF0cml4KH4uLCBkYXRhID0gY2JpbmQocHJzX21hdFssIGBQdF81ZS0wOGBdLCBjb3ZfbWF0KSkKICBwcnNfYW5ub3RhdGVkIDwtIGNwZy5hbm5vdGF0ZShkYXRhdHlwZSA9ICJhcnJheSIsIG9iamVjdCA9IG12YWx1ZXMsIGFuYWx5c2lzLnR5cGUgPSAiZGlmZmVyZW50aWFsIiwgZGVzaWduID0gZGVzaWduX3BycywgY29lZiA9IDIsIHdoYXQgPSAiTSIsIGFycmF5dHlwZSA9ICJFUElDIiwgZmRyID0gMC4wNSkKICBwcnNfZG1yX3JlcyA8LSBkbXJjYXRlKHByc19hbm5vdGF0ZWQsIGxhbWJkYSA9IDEwMDAsIEMgPSAyKQogIHJldHVybihTNF90b19kYXRhZnJhbWUocHJzX2Rtcl9yZXMpKQp9CmRtcl9jcm9zcyA8LSBydW5fZG1yY2F0ZShwcnNfbWF0LCBjb3ZfbWF0LCBtdmFsdWVzKQpkbXJfbWFsZXMgPC0gcnVuX2RtcmNhdGUocHJzX21hdF9tYWxlLGNvdl9tYXRfbWFsZSxtdmFsdWVzX21hbGUpCmRtcl9mZW1hbGUgPC0gcnVuX2RtcmNhdGUocHJzX21hdF9mZW1hbGUsY292X21hdF9mZW1hbGUsbXZhbHVlc19mZW1hbGUpCmBgYApgYGB7cn0KYW5ub3RhdGlvbiA8LSBtaW5maTo6Z2V0QW5ub3RhdGlvbihJbGx1bWluYUh1bWFuTWV0aHlsYXRpb25FUElDYW5uby5pbG0xMGI0LmhnMTkpCmBgYAojIyMgUGxvdHRpbmcgYWxsIERNUnMgey50YWJzZXQgLnRhYnNldC1mYWRlfQpgYGB7cn0KZ2V0X2Rtcl9lZmZlY3RzIDwtIGZ1bmN0aW9uKGRtciwgbGltbWFfcmVzLG12YWx1ZXMpIHsKICBkbXJfY29vcmQgPC0gc3RyX21hdGNoX2FsbChhcy5jaGFyYWN0ZXIoZG1yKSwgIihjaHIuKik6KFswLTldKiktKFswLTldKikiKVtbMV1dCiAgY3BncyA8LSBhcy5kYXRhLnRhYmxlKGFubm90YXRpb25bYW5ub3RhdGlvbiRjaHIgPT0gZG1yX2Nvb3JkWzJdICYgYW5ub3RhdGlvbiRwb3MgPj0gYXMubnVtZXJpYyhkbXJfY29vcmRbM10pICYgYW5ub3RhdGlvbiRwb3MgPD0gYXMubnVtZXJpYyhkbXJfY29vcmRbNF0pLCBdKQogIHJlcyA8LSBsaW1tYV9yZXNbY3Bncyxvbj1jKCJJRCI9Ik5hbWUiKSxub21hdGNoPTAgXQogIGRtcl8xIDwtIGFzLmRhdGEuZnJhbWUocmVzWyFpcy5uYShyZXMkbG9nRkMpLCBdKQogIGRtcl9tZXRoeSA8LSByZXNoYXBlMjo6bWVsdChsdW1pOjptMmJldGEobXZhbHVlc1tkbXJfMSRJRCwgXSksIHN0cmluZ3NBc0ZhY3RvcnMgPSBGQUxTRSkKICB0b19wbG90IDwtIG1lcmdlKGRtcl8xLCBkbXJfbWV0aHksIGJ5LnggPSAiSUQiLCBieS55ID0gIlZhcjEiKQogIHRvX3Bsb3QkRE1SIDwtIGRtcgogIHRvX3Bsb3QKfQpnZXRfZG1yX3JlcyA8LSBmdW5jdGlvbihkbXIsIGxpbW1hX3JlcykgewogIGRtcl9jb29yZCA8LSBzdHJfbWF0Y2hfYWxsKGFzLmNoYXJhY3RlcihkbXIpLCAiKGNoci4qKTooWzAtOV0qKS0oWzAtOV0qKSIpW1sxXV0KICBjcGdzIDwtIGFubm90YXRpb25bYW5ub3RhdGlvbiRjaHIgPT0gZG1yX2Nvb3JkWzJdICYgYW5ub3RhdGlvbiRwb3MgPj0gYXMubnVtZXJpYyhkbXJfY29vcmRbM10pICYgYW5ub3RhdGlvbiRwb3MgPD0gYXMubnVtZXJpYyhkbXJfY29vcmRbNF0pLCBdCiAgcmVzIDwtIGNiaW5kKGNwZ3MsIGxpbW1hX3Jlc1tjcGdzJE5hbWUsIF0pCiAgcmVzJERNUiA8LSBkbXIKICByZXR1cm4ocmVzKQp9CnBsb3RfZG1ycyA8LSBmdW5jdGlvbihkbXJfcmVzLGxpbW1hX3JlcyxwcnNfbWF0LCBtdmFsdWVzLCBjYXNlX2NvbnRyb2wpIHsKICB0b19wbG90IDwtIHJiaW5kbGlzdChsYXBwbHkoZG1yX3JlcyRjb29yZCwgZnVuY3Rpb24oZG1yKSBnZXRfZG1yX2VmZmVjdHMoZG1yLCBsaW1tYV9yZXMsbXZhbHVlcykpKSAlPiUKICAgIG11dGF0ZShTQ09SRTFfQVZHID0gc2NhbGUocHJzX21hdFttYXRjaChWYXIyLElEcyRJSUQpLGBQdF81ZS0wOGBdKSxQRCA9IFBEW21hdGNoKFZhcjIsSUlEKV0kUEQpICU+JQogICAgZmlsdGVyKCFpcy5uYShTQ09SRTFfQVZHKSkKICBmb3IoY3VyX2RtciBpbiB1bmlxdWUodG9fcGxvdCRETVIpKXsKICAgIGN1cl9wbG90IDwtIHRvX3Bsb3RbRE1SID09IGN1cl9kbXJdCiAgICBwMSA8LSBnZ3Bsb3QoY3VyX3Bsb3QgLCBhZXMocG9zLCB2YWx1ZSwgY29sb3IgPSBTQ09SRTFfQVZHKSkgKwogICAgICBnZW9tX3BvaW50KCkgKwogICAgICBzY2FsZV9jb2xvcl9ncmFkaWVudChsb3cgPSByZXYoY2FzZV9jb250cm9sKVsxXSwgaGlnaCA9IHJldihjYXNlX2NvbnRyb2wpWzJdKSArCiAgICAgIHRoZW1lX21pbmltYWwoKSArCiAgICAgIGxhYnModGl0bGUgPSB1bmlxdWUoY3VyX3Bsb3QkRE1SKSx4PSJQT1MiLHk9YnF1b3RlKCJNZXRoeWxhdGlvbiJ+YmV0YSksY29sb3IgPSAiTm9ybWFsaXplZCBQRCBQUlMiKSsKICAgICAgdGhlbWUoYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSA5MCkpCiAgICBwMiA8LSBnZ3Bsb3QoY3VyX3Bsb3QgJT4lIG11dGF0ZShQRCA9IGlmZWxzZShQRCA9PSAxLCAiQ0FTRSIsICJDT05UUk9MIikpLCBhZXMoZmFjdG9yKHBvcyksIHZhbHVlLCBjb2xvciA9IFBEKSkgKwogICAgICBnZW9tX2JveHBsb3QocG9zaXRpb24gPSBwb3NpdGlvbl9kb2RnZSgwLjc1KSkgKwogICAgICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzID0gY2FzZV9jb250cm9sKSArCiAgICAgIHRoZW1lX21pbmltYWwoKSArCiAgICAgIHN0YXRfc3VtbWFyeShhZXMoZ3JvdXAgPSBQRCksIGZ1biA9IG1lYW4sIGdlb20gPSAibGluZSIpICsKICAgICAgbGFicyh0aXRsZSA9IHVuaXF1ZShjdXJfcGxvdCRETVIpLHg9IlBPUyIseT1icXVvdGUoIk1ldGh5bGF0aW9uIn5iZXRhKSxjb2xvciA9ICJQRCBzdGF0dXMiKSArCiAgICAgIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gOTApKQogICAgcHJpbnQocDEpCiAgICBwcmludChwMikKICB9Cn0KCnBsb3RfZG1ycyhkbXJfY3Jvc3MsIHByc19hbm5vdCwgcHJzX21hdCxtdmFsdWVzLCByZXYoYygiZ3JheTgwIiwgImdyYXk0MCIpKSkKcGxvdF9kbXJzKGRtcl9tYWxlcyxwcnNfYW5ub3RfbWFsZSxwcnNfbWF0X21hbGUsIG12YWx1ZXNfbWFsZSxyZXYoYygibGlnaHQgYmx1ZSIsICJsaWdodGJsdWU0IikpKQpwbG90X2RtcnMoZG1yX2ZlbWFsZSwgcHJzX2Fubm90X2ZlbWFsZSwgcHJzX21hdF9mZW1hbGUsIG12YWx1ZXNfZmVtYWxlLHJldihjKCJsaWdodHBpbmsiLCAibGlnaHRwaW5rNCIpKSkKYGBgCgoKYGBge3J9CmRtcl9jcm9zcwpkbXJfbWFsZXMKZG1yX2ZlbWFsZQpzYXZlKGxpc3Q9YygiZG1yX2Nyb3NzIiwiZG1yX21hbGVzIiwiZG1yX2ZlbWFsZSIpLGZpbGU9InByc19kbXJfbmFsbHMuUkRhdGEiKQpgYGAK